%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes diagonal blocks of the coefficient matrix.
% Note that this should be used ONLY for an index set ind that is
% entirely contained within one contour.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = get_A_single_diag(C,ind,kh)

ntot     = size(C,2);
nloc     = length(ind);

% Y_g1(i,j) = x(j) , X_g1(i,j) = x(i)
[Y_g1,   X_g1   ] = meshgrid(C(1,ind), C(1,ind));
% Y_g2(i,j) = y(j) , X_g2(i,j) = y(i)
[Y_g2,   X_g2   ] = meshgrid(C(4,ind), C(4,ind));
% Y_dg1(i,j) = dx(j) , X_dg1(i,j) = dx(i)
[Y_dg1,  X_dg1  ] = meshgrid(C(2,ind), C(2,ind));
% Y_dg2(i,j) = dy(j) , X_dg2(i,j) = dy(i)
[Y_dg2,  X_dg2  ] = meshgrid(C(5,ind), C(5,ind));

ima = sqrt(-1);
% tangent vector : t = (dx,dy)
% -> normal vector n = (dy,-dx) / |t| 
% nn1(i,j) = nx(j)
nn1 = ( Y_dg2./sqrt(Y_dg1.*Y_dg1 + Y_dg2.*Y_dg2));
% nn2(i,j) = ny(j)
nn2 = (-Y_dg1./sqrt(Y_dg1.*Y_dg1 + Y_dg2.*Y_dg2));
% dd(i,j) = delta(i,j) + |x(j) - x(i)|
dd  = sqrt((Y_g1 - X_g1).^2 + (Y_g2 - X_g2).^2) + eye(nloc);

K   = funcK(kh,nn1,nn2,(Y_g1 - X_g1),(Y_g2 - X_g2),dd);

% for i !=j
% A(i,j) = h |n(j)| K(x(j) - x(i))
% for i = j
% A(i,j) = -2i

A   = K.*sqrt(Y_dg1.^2 + Y_dg2.^2);

MU6 = [ 0.2051970990601250e1 + 0.2915391987686505e1;...
        -0.7407035584542865e1 - 0.8797979464048396e1;...
        0.1219590847580216e2 + 0.1365562914252423e2;...
        -0.1064623987147282e2 - 0.1157975479644601e2;...
        0.4799117710681772e1 + 0.5130987287355766e1;...
        -0.8837770983721025   - 0.9342187797694916];

% distance between two indices
% D(i,j) = |ind(i) - ind(j)|
D = abs(ind'*ones(1,length(ind)) - ones(length(ind),1)*ind);
% take into account periodicity 
D = abs(mod(D + 6, ntot) - 6);
% extract pairs of indices that are closer than 6 discretization points 
[ii1, ii2] = find(D.*(D<7));
for j = 1:length(ii1)
    i1 = ii1(j);
    i2 = ii2(j);
    d  = D(i1,i2);
    A(i1,i2) = (1 + MU6(d))*A(i1,i2);
end

% diagonal terms
for j = 1:nloc
    A(j,j) = - 2*ima;
end

return
